!// Dmitry.Konovalov@jcu.edu.au Dmitry.A.Konovalov@gmail.com, 2013

module vec_mod; 
use dbg_mod; 
implicit none;
!private;
logical, save :: VEC_TRACE = .FALSE.;
integer, parameter :: DUFF_SIZE = 8;
integer, parameter :: DUFF_SIZE_16 = 16;

contains; 




subroutine vec_dbg(info, v, trace);   
  character(len=1024) :: mssg;
  character(LEN=*), intent(IN) :: info;   
  real(DP), intent(IN) :: v(:); 
  logical, optional, intent(in) :: trace;  
  integer :: i, n, mini, maxi;
  if (present(trace) .AND. (.NOT. trace)) RETURN;
  
  call vec_range(mini, maxi, v, n);
  write(mssg,*) info, "(", TRIM(int_to_str(mini)), &
    & ":", TRIM(int_to_str(maxi)), ")="; 
  if (n <= VEC_DBG_VIEW_N) then
    write(*,*) TRIM(mssg),  (v(i), i = mini, maxi);
  else
    n = VEC_DBG_VIEW_N / 2;
    write(*,*) TRIM(mssg), &
    & (v(i), i = mini, mini + n - 1), ", ..., ", &
    & (v(i), i = maxi - n + 1, maxi);
  end if
end subroutine


subroutine veci_dbg(info, iv, trace);   
  character(len=1024) :: mssg;
  character(LEN=*), intent(IN) :: info;   
  integer, intent(IN) :: iv(:); 
  logical, optional, intent(in) :: trace;  
  integer :: i, n, mini, maxi;
  if (present(trace) .AND. (.NOT. trace)) RETURN;
  
  call veci_range(mini, maxi, iv, n);
  write(mssg,*) info, "(", TRIM(int_to_str(mini)), &
    & ":", TRIM(int_to_str(maxi)), ")="; 
  if (n <= VEC_DBG_VIEW_N) then
    write(*,*) TRIM(mssg),  (iv(i), i = mini, maxi);
  else
    n = VEC_DBG_VIEW_N / 2;
    write(*,*) TRIM(mssg), &
    & (iv(i), i = mini, mini + n - 1), ", ..., ", &
    & (iv(i), i = maxi - n + 1, maxi);
  end if
end subroutine



SUBROUTINE vec_init_res(res, x);    
  real(DP), intent(OUT) :: res(:);      ! vector-x, vector-y
  real(DP), intent(IN)  :: x(:);      
  integer :: mini, maxi, minr, maxr;
  res = 0; ! cleanup first. size of res may be > than of vx
  call vec_range(mini, maxi, x); 
  call vec_range(minr, maxr, res); 
  call dbg_assert_true("minr == mini", minr == mini);
  call dbg_assert_true("maxi <= maxr", maxi <= maxr);
end SUBROUTINE

SUBROUTINE vec_assert_math(v, v2);   
  real(DP), intent(IN)  :: v(:), v2(:);      
  integer :: mini, maxi, min2, max2;
  call vec_range(mini, maxi, v); 
  call vec_range(min2, max2, v2); 
  call dbg_assert_true("mini == min2", mini == min2);
  call dbg_assert_true("maxi == max2", maxi == max2);
end SUBROUTINE

subroutine vec_range(mini, maxi, v, nv);  
  real(DP), intent(IN) :: v(:); 
  integer, intent(OUT) :: mini, maxi;
  integer, intent(OUT), optional :: nv;
  mini = LBOUND(v, 1);
  maxi = UBOUND(v, 1);
  if (present(nv)) nv = SIZE(v, 1);
end SUBROUTINE

subroutine veci_range(mini, maxi, iv, nv);  
  integer, intent(IN) :: iv(:); 
  integer, intent(OUT) :: mini, maxi;
  integer, intent(OUT), optional :: nv;
  mini = LBOUND(iv, 1);
  maxi = UBOUND(iv, 1);
  if (present(nv)) nv = SIZE(iv, 1);
end SUBROUTINE

subroutine vec_div_SAFE(res, top, v);   
  real(DP), intent(OUT) :: res(:);  
  real(DP), intent(IN) :: top, v(:);      
  integer :: minr, maxr, mini, maxi, i, ir;
  call vec_range(minr, maxr, res);
  call vec_range(mini, maxi, v);
  res = D0;
  ir = minr;
  do i = mini, maxi
    if (ir > maxr  .OR.  i > maxi)  return;   
    if (dbg_is_zero(v(i))) then
      res(ir) = D0;
    else
      res(ir) = top / v(i);
    endif
    ir = ir + 1;
  enddo  
end subroutine

subroutine vec_copy_SAFE(res, v);   
  real(DP), intent(OUT) :: res(:);  
  real(DP), intent(IN) :: v(:);      
  integer :: minr, maxr, mini, maxi, i, ir;
  call vec_range(minr, maxr, res);
  call vec_range(mini, maxi, v);
  res = D0;
  ir = minr;
  do i = mini, maxi
    if (ir > maxr  .OR.  i > maxi)  return;
    res(ir) = v(i);
    ir = ir + 1;
  enddo  
end subroutine

subroutine veci_copy_SAFE(res, iv);   
  integer, intent(OUT) :: res(:);  
  integer, intent(IN) :: iv(:);      
  integer :: minr, maxr, mini, maxi, i, ir;
  call dbg_str("-->veci_copy_SAFE()", VEC_TRACE);
  call veci_range(minr, maxr, res);
  call veci_range(mini, maxi, iv);
  res = 0;
  ir = minr;
  do i = mini, maxi
    if (ir > maxr  .OR.  i > maxi)  return;
    res(ir) = iv(i);
    ir = ir + 1;
  enddo  
end subroutine



real(DP) function vec_dot_prod(v, v2) result(res);  implicit none;  
  real(DP), dimension(:), intent(in) :: v, v2;  
  call dbg_str("-->vec_dot_prod()", VEC_TRACE);
  res = DOT_PRODUCT(v, v2);
  return
end function

real(DP) function vec_dot_prod2(v, v2) result(res);  implicit none;  
  real(DP), dimension(:), intent(in) :: v, v2;  
!  call dbg_str("-->vec_dot_prod2()", VEC_TRACE);
  res = SUM(v * v2);
  return
end function

real(DP) function vec_dot_prod3(v, v2, v3) result(res);  implicit none;  
  real(DP), dimension(:), intent(in) :: v, v2, v3;  
!  call dbg_str("-->vec_dot_prod3()", VEC_TRACE);
  res = SUM(v * v2 * v3);
  return
end function

real(DP) function vec_dot_prod4(v, v2, v3, v4) result(res);  implicit none;  
  real(DP), dimension(:), intent(in) :: v, v2, v3, v4;  
 ! call dbg_str("-->vec_dot_prod4()", VEC_TRACE);
  res = SUM(v * v2 * v3 * v4);
  return
end function



! assert two vectors are the same
subroutine vec_assert_abs_diff_LE(info, v, v2, maxerr);  
  character(len=1024) :: mssg;
  character(LEN=*), intent(IN) :: info;
  real(DP), intent(IN) :: v(:), v2(:), maxerr;
  integer :: i, maxi, mini, min2, max2;
  character(LEN=64) :: name = "vec_assert_abs_diff_LE";
  call dbg_in(name);
  call vec_range(mini, maxi, v);
  call vec_range(min2, max2, v2);
  call dbg_assert_true("mini == min2", mini == min2);
  call dbg_assert_true("maxi == max2", maxi == max2);
  do i = mini, maxi
    write (mssg, *) "-->vec_assert_abs_diff_LE(info=", info, ", i = ", i;
!    call dbg_dp(TRIM(mssg), v(i));
    !call dbg_int("i=", i);
    !call dbg_dp("   v(i)=", v(i));
    !call dbg_dp("  v2(i)=", v2(i));
    !call dbg_dp("abs err=", ABS(v(i) - v2(i)));
    call dbg_assert_abs_diff_LE(TRIM(mssg), v(i), v2(i), maxerr);
  end do
end subroutine



!================================ DO NOT USE BELOW ============




real(DP) function vec_dot_do(v, v2) result(res);  implicit none;  
  real(DP), dimension(:), intent(in) :: v, v2;   
  integer :: i;
  call dbg_str("-->vec_dot_do()", VEC_TRACE);
  res = 0.;
  do i = LBOUND(v, 1), UBOUND(v, 1)
    res = res + v(i) * v2(i);
  end do 
  return
end function

real(DP) function vec_dot_do3(v, v2, v3) result(res);  implicit none;  
  real(DP), dimension(:), intent(in) :: v, v2, v3;   
  integer :: i;
  call dbg_str("-->vec_dot_do3()", VEC_TRACE);
  res = 0.;
  do i = LBOUND(v, 1), UBOUND(v, 1)
    res = res + v(i) * v2(i) * v3(i);
  end do 
  return
end function


real(DP) function vec_dot_duff(v, v2) result(res);  implicit none;  
  real(DP), dimension(:), intent(in) :: v, v2;   
  integer :: i, m, i8, min_i, max_i, n, start_m, start_i8;
  call dbg_str("-->vec_dot_duff()", VEC_TRACE);
  min_i = LBOUND(v, 1); max_i = UBOUND(v, 1);
  res = 0.;
  i = min_i;
  n = max_i - min_i + 1;         !call dbg_int("n=", n);
  start_m = MOD(n, DUFF_SIZE);   !call dbg_int("start_m=", start_m);
  do m = start_m, 1, -1 
    res = res + v(i) * v2(i);  i = i + 1;      
    !call dbg_int("i=", i);
  end do
  start_i8 = n / DUFF_SIZE;
  !call dbg_int("start_i8=", start_i8);
  !call dbg_int("before do i=", i);
  do i8 = start_i8, 1, -1
    res = res + v(i) * v2(i);      i = i + 1; !1
    res = res + v(i) * v2(i);      i = i + 1; !2
    res = res + v(i) * v2(i);      i = i + 1; !3
    res = res + v(i) * v2(i);      i = i + 1; !4
    res = res + v(i) * v2(i);      i = i + 1; !5
    res = res + v(i) * v2(i);      i = i + 1; !6
    res = res + v(i) * v2(i);      i = i + 1; !7
    res = res + v(i) * v2(i);      i = i + 1; !8
  end do
  !call dbg_int("after do i=", i);
  return
end function

real(DP) function vec_dot_duff3(v, v2, v3) result(res);  implicit none;  
  real(DP), dimension(:), intent(in) :: v, v2, v3;   
  integer :: i, m, i8, min_i, max_i, n, start_m, start_i8;
  call dbg_str("-->vec_dot_duff3()", VEC_TRACE);
  min_i = LBOUND(v, 1); max_i = UBOUND(v, 1);
  res = 0.;
  i = min_i;
  n = max_i - min_i + 1;         !call dbg_int("n=", n);
  start_m = MOD(n, DUFF_SIZE);   !call dbg_int("start_m=", start_m);
  do m = start_m, 1, -1 
    res = res + v(i) * v2(i)* v3(i);  i = i + 1;      
    !call dbg_int("i=", i);
  end do
  start_i8 = n / DUFF_SIZE;
  !call dbg_int("start_i8=", start_i8);
  !call dbg_int("before do i=", i);
  do i8 = start_i8, 1, -1
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !1
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !2
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !3
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !4
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !5
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !6
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !7
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !8
  end do
  !call dbg_int("after do i=", i);
  return
end function

real(DP) function vec_dot_duff3c(v, v2, v3) result(res);  implicit none;  
  real(DP), dimension(:), intent(in) :: v, v2, v3;   
  integer :: i, m, i8, min_i, max_i, n, start_m, start_i8;
  call dbg_str("-->vec_dot_duff3c()", VEC_TRACE);
  min_i = LBOUND(v, 1); max_i = UBOUND(v, 1);
  res = 0.;
  i = min_i;
  n = max_i - min_i + 1;         !call dbg_int("n=", n);
  start_m = MOD(n, DUFF_SIZE_16);   !call dbg_int("start_m=", start_m);
  do m = start_m, 1, -1 
    res = res + v(i) * v2(i)* v3(i);  i = i + 1;      
    !call dbg_int("i=", i);
  end do
  start_i8 = n / DUFF_SIZE_16;
  !call dbg_int("start_i8=", start_i8);
  !call dbg_int("before do i=", i);
  do i8 = start_i8, 1, -1
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !1
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !2
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !3
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !4
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !5
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !6
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !7
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !8
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !9
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !10
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !11
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !12
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !13
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !14
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !15
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !16
  end do
  !call dbg_int("after do i=", i);
  return
end function

real(DP) function vec_dot_duff3b(v, v2, v3) result(res);  implicit none;  
  real(DP), dimension(:), intent(in) :: v, v2, v3;   
  integer :: i, m, i8, min_i, max_i, n, start_m, start_i8;
  call dbg_str("-->vec_dot_duff3b()", VEC_TRACE);
  min_i = LBOUND(v, 1); max_i = UBOUND(v, 1);
  res = 0.;
  i = min_i;
  n = max_i - min_i + 1;         !call dbg_int("n=", n);
  start_m = MOD(n, DUFF_SIZE_16);   !call dbg_int("start_m=", start_m);
  do m = start_m, 1, -1 
    res = res + v(i) * v2(i)* v3(i);  i = i + 1;      
    !call dbg_int("i=", i);
  end do
  start_i8 = n / DUFF_SIZE_16;
  !call dbg_int("start_i8=", start_i8);
  !call dbg_int("before do i=", i);
  do i8 = start_i8, 1, -1
    res = res + v(i) * v2(i) * v3(i);      !i = i + 1; !1
    res = res + v(i+1) * v2(i+1) * v3(i+1);      !i = i + 1; !2
    res = res + v(i+2) * v2(i+2) * v3(i+2);      !i = i + 1; !3
    res = res + v(i+3) * v2(i+3) * v3(i+3);      !i = i + 1; !4
    res = res + v(i+4) * v2(i+4) * v3(i+4);      !i = i + 1; !5
    res = res + v(i+5) * v2(i+5) * v3(i+5);      !i = i + 1; !6
    res = res + v(i+6) * v2(i+6) * v3(i+6);      !i = i + 1; !7
    res = res + v(i+7) * v2(i+7) * v3(i+7);      
    res = res + v(i+8) * v2(i+8) * v3(i+8);      
    res = res + v(i+9) * v2(i+9) * v3(i+9);      
    res = res + v(i+10) * v2(i+10) * v3(i+10);      
    res = res + v(i+11) * v2(i+11) * v3(i+11);      
    res = res + v(i+12) * v2(i+12) * v3(i+12);      
    res = res + v(i+13) * v2(i+13) * v3(i+13);      
    res = res + v(i+14) * v2(i+14) * v3(i+14);      
    res = res + v(i+15) * v2(i+15) * v3(i+15);      
    i = i + DUFF_SIZE_16; 
  end do
  !call dbg_int("after do i=", i);
  return
end function
end module 
